/*
 * CUDA_vec.cu
 *
 *  Created on: Mar 22, 2011
 *      Author: curt
 */
#include "CUDA_vec.h"
#include "Operations.h"
#include "Reduce.h"
#include "utility.h"
#include <stdio.h>

const static int BLOCK_SIZE = 256;

template<typename T>
CUDA_vec<T>::CUDA_vec(){
	size = 0;
	data = NULL;
}

template <typename T>
void CUDA_vec<T>::set_constant(T constant)
{
	if(size!=0)
	{
		T * array = new T[size];
		for(int i=0; i<size; i++)
			array[i] = constant;
		cudaMemcpy(data, array, sizeof(T) * size, cudaMemcpyHostToDevice);
		delete [] array;
	}


}
template <typename T>
void CUDA_vec<T>::extractFromDevice (T * host_data) const
{
	cudaMemcpy(host_data, data,  sizeof(T) * size, cudaMemcpyDeviceToHost);
}
template <typename T>
void CUDA_vec<T>::copyToDevice(T * host_data, int host_size){
	if(host_size <= size)
	{
		cudaMemcpy(data, host_data,  sizeof(T) * host_size, cudaMemcpyHostToDevice);
	}
	else
	{
		printf("Error: trying to copy more data to device vector than allocated\n");
		exit(1);
	}
}
template <typename T>
void CUDA_vec<T>::destroy(){
	if(data!=NULL)
	{
		cudaFree(data);
		data = NULL;
	}
	else
	{
		printf("Tried to free a null vector\n");
	}
}
template <typename T>
CUDA_vec<T>::CUDA_vec(int s) {
	size = s;
	cudaMalloc(&data, sizeof(T) * size);
	print_CUDA_err_msg("Vector initialization");
	set_constant(0);
}
template <typename T>
CUDA_vec<T>::CUDA_vec(const CUDA_vec<T> & b)
{
	size = b.getSize();
	cudaMalloc(&data, sizeof(T) * size);
	print_CUDA_err_msg("Copy constructor: Before memcopy");
	cudaMemcpy(data, b.raw_pointer(), sizeof(T) * size, cudaMemcpyDeviceToDevice);
	print_CUDA_err_msg("Copy constructor: After memcopy");
}

template <typename T>
const T * CUDA_vec<T>::raw_pointer() const { return data; }

template <typename T>
CUDA_vec<T>::~CUDA_vec() {
	destroy();
}
template<typename T>
T * CUDA_vec<T>::raw_pointer(){ return data; }

template<typename T>
void CUDA_vec<T>::copyNewArray(const T * array, int s)
{
	if(data!=NULL)
		cudaFree(data);
	size = s;
	cudaMalloc(&data, sizeof(T) * size);
	copySameArray(array, size);
}
template<typename T>
void CUDA_vec<T>::copySameArray(const T * array, int s)
{
	if(data!=NULL && array!=NULL && s<= size)
		cudaMemcpy(data, array, sizeof(T) * s, cudaMemcpyDeviceToDevice);
	else
	{
		printf("data %d array %d s:%d size:%d \n", data, array, s, size);
		printf("Error: either null arrays or invalid copy sizes\n");
		exit(1);
	}
}

template<typename T>
CUDA_vec<T> & CUDA_vec<T>::operator=(const CUDA_vec<T> & other)
{
	if(&other!=this)
	{
		data = other.data;
		size = other.size;
	}
	return *this;
}
template<typename T>
void CUDA_vec<T>::copyValues(const CUDA_vec<T> & other)
{
	print_CUDA_err_msg("Before memory allocation");
	copySameArray(other.raw_pointer(), other.getSize());
	print_CUDA_err_msg("CUDA_vec memory allocation");
}

template<typename T>
__global__ void axpby(T a, T b, T * x, const T * y, int length)
{
	int i = blockIdx.x * blockDim.x + threadIdx.x;
	if(i <  length)
		x[i] = a*x[i] + b*y[i];
}

template<typename T>
__global__ void recip(T * x, int length)
{
	int i = blockIdx.x * blockDim.x + threadIdx.x;
	if(i <  length)
		x[i] = 1/x[i];
}

template<typename T>
__global__ void ptwise_mult(T * x, const T * y, int length)
{
	int i = blockIdx.x * blockDim.x + threadIdx.x;
	if(i <  length)
		x[i] = x[i] * y[i];
}

template<typename T>
CUDA_vec<T> & CUDA_vec<T>::AXPBY(T a, T b, const CUDA_vec<T> & y)
{
	dim3 blockSize(BLOCK_SIZE);
	dim3 gridSize( size / BLOCK_SIZE + 1 );
	axpby<T><<<gridSize,blockSize>>>(a, b, data, y.raw_pointer(), size);
	return *this;
}

template<typename T>
CUDA_vec<T> & CUDA_vec<T>::reciprocal()
{
	dim3 blockSize(BLOCK_SIZE);
	dim3 gridSize( size / BLOCK_SIZE + 1 );
	recip<T><<<gridSize,blockSize>>>(data, size);
	return *this;
}
template<typename T>
CUDA_vec<T> & CUDA_vec<T>::operator*=(const CUDA_vec<T> & other)
{
	dim3 blockSize(BLOCK_SIZE);
	dim3 gridSize( size / BLOCK_SIZE + 1 );
	ptwise_mult<T><<<gridSize,blockSize>>>(data, other.raw_pointer(), size);
	print_CUDA_err_msg("After *=");
	return *this;
}
template<typename T>
T CUDA_vec<T>::dotProd(const CUDA_vec<T> & other) const
{
	print_CUDA_err_msg("Before temp allocation");
	CUDA_vec<T> temp(*this);
	print_CUDA_err_msg("Before *=");
	temp *= other;
	plus<T> op;
	T dotprod = CUDA_reduce(temp, op);
	return dotprod;

}

template class CUDA_vec<float>;
//template class CUDA_vec<int>;
//template class CUDA_vec<unsigned>;
//template class CUDA_vec<unsigned long>;
//template class CUDA_vec<long>;
